package mat;

import java.io.BufferedOutputStream;
import java.io.DataOutputStream;
import java.io.IOException;
import java.io.OutputStream;
import java.io.PrintStream;
import java.math.BigDecimal;
import java.math.MathContext;
import java.math.RoundingMode;
import java.nio.ByteBuffer;
import java.nio.DoubleBuffer;
import java.nio.LongBuffer;
import java.util.ArrayList;
import java.util.Arrays;
import java.util.Comparator;
import java.util.Formatter;
import java.util.LinkedList;
import java.util.List;
import java.util.concurrent.ExecutorService;
import java.util.concurrent.Future;
import java.util.concurrent.LinkedBlockingQueue;
import java.util.concurrent.ThreadPoolExecutor;
import java.util.concurrent.TimeUnit;
import java.util.concurrent.atomic.AtomicInteger;

import org.apache.log4j.Logger;

import Jama.Matrix;
import Jama.QRDecomposition;


public class Mat
{
	public static Logger log = Logger.getLogger(Mat.class);
	public static boolean verbose = false;
	public static ThreadPoolExecutor executorService;
	
	/*
	 * Number of threads for "simple" parallelized methods
	 * FIXME
	 * Split Mat into a SingleThreadedMat and MultiThreadedMat.
	 * The value of simpleMethodMaxThreads shouldn'be set empirically,
	 * but derived from e.g. the size of the thread pool. Note
	 * that if Mat instances run concurrently in several threads, they
	 * will use the same (i.e. static) Mat thread pool for multiply()
	 * and parallelized functions. For simple functions, it makes no
	 * sense to eat up the entire pool per thread since this would
	 * block other threads from using the pool. Long story short: 
	 * simpleMethodMaxThreads should be significantly smaller than
	 * the size of the thread pool if run in a multithreaded environment -
	 * which is the case for scengen.
	 */
	public static int simpleMethodMaxThreads = 4;
	public static boolean testScaling = false;

	private String name;
	private final double[][] data;
	private final int rowDimension;
	private final int colDimension;
	
	private boolean isInner;
	
	public enum FileFormat { DEFAULT, REVERSE_BYTE_ORDER }
	private FileFormat outputFileFormat = FileFormat.DEFAULT;

	public static int ROWS = 2;
	public static int COLS = 1;
	
	private static AtomicInteger asyncWriteInProgress = new AtomicInteger();
	
	public void setOutputFileFormat( FileFormat off )
	{
		outputFileFormat = off;
	}
	
	public static ExecutorService setupThreadPool(int maxThreads)
	{
		executorService = new ThreadPoolExecutor(maxThreads, maxThreads,
                0l, TimeUnit.MILLISECONDS,
                new LinkedBlockingQueue<Runnable>(),
                new MatThreadFactory());			
        return executorService;
	}
	
	public static int getMaxThreads()
	{
		return executorService == null ? 0: executorService.getMaximumPoolSize();	
	}
	
	public static int getAsyncWriteInProgress()
	{
		return asyncWriteInProgress.get();
	}
	
	public static boolean isAsyncWriteInProgress()
	{
		return asyncWriteInProgress.get() > 0;
	}
	
	public static int wait4AsyncWriteInProgress(long timeoutInSeconds) throws InterruptedException
	{
		long t = System.currentTimeMillis() + timeoutInSeconds*1000L;
		while (System.currentTimeMillis() < t)
		{
			if (getAsyncWriteInProgress() < 1)
				return 0;
			
			Thread.sleep(5000);
		}
		
		int threads = getAsyncWriteInProgress();
		log.warn("Still " + threads + "async Mat writers running!");
		return threads;
	}
	

	public static String concat(Object... objs)
	{
		StringBuilder sb = new StringBuilder();
		for (Object obj : objs)
			sb.append(obj);
		return sb.toString();
	}

	public static Mat readArray(String name, String fileName, int rows, int cols) throws IOException, InterruptedException
	{
		BinaryMatrixReader reader = new BinaryMatrixReader(name, fileName, rows, cols);
		reader.process();
		Mat m = reader.getMatrix();

		if (verbose)
			m.print();

		return m;
	}

	public static Mat load(String name, String fileName) throws IOException, FileParserException
	{
		TextMatrixReader reader = new TextMatrixReader(name, fileName);
		reader.process();
		Mat m = reader.getMatrix();

		if (verbose)
			m.print();

		return m;
	}

	public static Mat diag(double[] diag)
	{
		double[][] data = new double[diag.length][diag.length];
		for (int i = 0; i < data.length; i++)
		{
			data[i][i] = diag[i];
		}
		return new Mat(data);
	}
	
	public static Mat diag(String name, int dim, double diag)
	{
		double[][] data = new double[dim][dim];
		for (int i = 0; i < data.length; i++)
		{
			data[i][i] = diag;
		}
		return new Mat(name, data);
	}
	
	public static Mat diag(String name, int rows, int cols, double diag)
	{
		double[][] data = new double[rows][cols];
		for (int i = 0; i < rows; i++)
		{
			for (int j = 0; j < cols; j++)
			{
				if (i == j)
					data[i][j] = diag;
			}
		}
		return new Mat(name, data);
	}


	public static Mat ones(int rows, int cols)
	{
		Mat ones = fillNmatix(rows, cols, 1.0);
		ones.setName("ones(" + rows + "," + cols + ")");
		return ones;			
	}
	
	
	public Mat compareEQ(Mat m)
	{
		double[][] result = createResult();
		double[][] mData = m.getData();
		for (int i = 0; i < rowDimension; i++)
		{
			double[] a = data[i];
			double[] b = mData[i];
			double[] c = result[i];

			for (int j = 0; j < colDimension; j++)
			{
				c[j] = (b[j] == a[j]) ? 1 : 0;
			}
		}
		
		return innerOrNew("(" + this.name + "==" + m.name + ")", result);
	}
	public boolean equals(Mat m)
	{
		double[][] mData = m.getData();
		if(rowDimension!=m.rowDimension||colDimension!=m.colDimension)
			return false;
		for (int i = 0; i < rowDimension; i++)
		{
			double[] a = data[i];
			double[] b = mData[i];
			for (int j = 0; j < colDimension; j++)
			{
				if(b[j] != a[j])
					return false;
			}
		}
		return true;
	}
	
	// v > m
	public Mat compareGT(double v)
	{
		double[][] result = createResult();
		for (int i = 0; i < rowDimension; i++)
		{
			double[] a = data[i];
			double[] c = result[i];
			for (int j = 0; j < colDimension; j++)
			{
				c[j] = v > a[j] ? 1.0 : 0.0;
			}
		}
		
		return innerOrNew("(" + v + ">" + name + ")", result);
	}
	
	public Mat absFloor(double threshold, double floor)
	{
		double[][] result = createResult();
		for (int i = 0; i < rowDimension; i++)
		{
			double[] a = data[i];
			double[] c = result[i];
			for (int j = 0; j < colDimension; j++)
			{
				double d = a[j];
				if (Math.abs(d) < threshold)
					d = floor;
				
				c[j] = d;
			}
		}

		return innerOrNew("absFloor(" + name + ")", result);
	}


	// min(v,m)
	public Mat min(double v)
	{
		double[][] result = createResult();
		for (int i = 0; i < rowDimension; i++)
		{
			double[] a = data[i];
			double[] c = result[i];
			for (int j = 0; j < colDimension; j++)
			{
				c[j] = Math.min(v, a[j]);
			}
		}

		return innerOrNew("min(" + name + ")", result);
	}

	// max(v,m)
	public Mat max(double v)
	{
		double[][] result = createResult();
		for (int i = 0; i < rowDimension; i++)
		{
			double[] a = data[i];
			double[] c = result[i];
			for (int j = 0; j < colDimension; j++)
			{
				c[j] = Math.max(v, a[j]);
			}
		}

		return innerOrNew("min(" + name + ")", result);
	}

	// x<t, x:=mat, t scalar
	public Mat lt(double v)
	{
		double[][] result = createResult();
		for (int i = 0; i < rowDimension; i++)
		{
			double[] a = data[i];
			double[] c = result[i];
			for (int j = 0; j < colDimension; j++)
			{
				c[j] = a[j] < v ? 1.0 : 0.0;
			}
		}

		return innerOrNew(name + "<" + v, result);
	}

	// x>=t, x:=mat, t scalar
	public Mat gteq(double v)
	{
		double[][] result = createResult();
		for (int i = 0; i < rowDimension; i++)
		{
			double[] a = data[i];
			double[] c = result[i];
			for (int j = 0; j < colDimension; j++)
			{
				c[j] = a[j] >= v ? 1.0 : 0.0;
			}
		}

		return innerOrNew(name + ">=" + v, result);
	}

	// x>t, x:=mat, t scalar
	public Mat gt(double v)
	{
		double[][] result = createResult();
		for (int i = 0; i < rowDimension; i++)
		{
			double[] a = data[i];
			double[] c = result[i];
			for (int j = 0; j < colDimension; j++)
			{
				c[j] = a[j] > v ? 1.0 : 0.0;
			}
		}

		return innerOrNew(name + ">=" + v, result);
	}

	// x==t, x:=mat, t scalar
	public Mat eq(double v)
	{
		double[][] result = createResult();
		for (int i = 0; i < rowDimension; i++)
		{
			double[] a = data[i];
			double[] c = result[i];
			for (int j = 0; j < colDimension; j++)
			{
				c[j] = a[j] == v ? 1.0 : 0.0;
			}
		}

		return innerOrNew(name + "==" + v, result);
	}

	public double min()
	{
		double min = Double.MAX_VALUE;

		for (int i = 0; i < rowDimension; i++)
		{
			double[] row = data[i];
			for (int j = 0; j < colDimension; j++)
			{
				double d = row[j];
				if (d < min)
					min = d;
			}
		}

		return min;
	}

	public double max()
	{
		double max = Double.MIN_VALUE;

		for (int i = 0; i < rowDimension; i++)
		{
			double[] row = data[i];
			for (int j = 0; j < colDimension; j++)
			{
				double d = row[j];
				if (d > max)
					max = d;
			}
		}

		return max;
	}
	
	public int maxColumnIndex()
	{
		double max = Double.MIN_VALUE;
		int col=-1;
		for (int i = 0; i < rowDimension; i++)
		{
			double[] row = data[i];
			for (int j = 0; j < colDimension; j++)
			{
				double d = row[j];
				if (d > max)
				{
					max = d;
					col=j;
				}
			}
		}
		return col;
	}


	public Mat pmax(Mat m)
	{
		double[][] result = createResult();
		double[][] mdata = m.getData();

		for (int i = 0; i < rowDimension; i++)
		{
			double[] a = data[i];
			double[] b = mdata[i];
			double[] c = result[i];
			for (int j = 0; j < colDimension; j++)
			{
				c[j] = Math.max(a[j], b[j]);
			}
		}

		return innerOrNew("pmax(" + name + ", " + m.name + ")", result);
	}

	public Mat pmin(Mat m)
	{
		double[][] result = createResult();
		double[][] mdata = m.getData();

		for (int i = 0; i < rowDimension; i++)
		{
			double[] a = data[i];
			double[] b = mdata[i];
			double[] c = result[i];
			for (int j = 0; j < colDimension; j++)
			{
				c[j] = Math.min(a[j], b[j]);
			}
		}

		return innerOrNew("pmin(" + name + ", " + m.name + ")", result);
	}

	public Mat pmin(double d)
	{
		double[][] result = createResult();

		for (int i = 0; i < rowDimension; i++)
		{
			double[] a = data[i];
			double[] c = result[i];
			for (int j = 0; j < colDimension; j++)
			{
				c[j] = Math.min(a[j], d);
			}
		}

		return innerOrNew("pmin(" + name + ", " + d + ")", result);
	}

	public Mat pmax(double d)
	{
		double[][] result = createResult();

		for (int i = 0; i < rowDimension; i++)
		{
			double[] a = data[i];
			double[] c = result[i];
			for (int j = 0; j < colDimension; j++)
			{
				c[j] = Math.max(a[j], d);
			}
		}

		return innerOrNew("pmax(" + name + ", " + d + ")", result);
	}

	public Mat pow(double v)
	{
		double[][] result = createResult();
			
		MatRowWorker job = new MatRowWorker("(" + v + ".^" + name + ")", result, data, v) {

			public void run()
			{
				for (int i = iStart; i < iEnd; i++)
				{
					double[] a = A[i];
					double[] c = result[i];
					for (int j = 0; j < a.length; j++)
					{
						c[j] = Math.pow(value, a[j]);
					}
				}			
			}
		};
		return executeMatRowWorkers(job, simpleMethodMaxThreads);
	}

	public double sum()
	{
		double sum = 0.;
		for (int i = 0; i < rowDimension; i++)
		{
			double[] a = data[i];
			for (int j = 0; j < colDimension; j++)
				sum += a[j];
		}

		return sum;
	}

	public Mat sumColumnWise()
	{
		double[][] sum = new double[1][colDimension];

		for (int i = 0; i < rowDimension; i++)
		{
			double[] a = data[i];
			double[] c = sum[0];
			for (int j = 0; j < colDimension; j++)
				c[j] += a[j];
		}

		return new Mat("sumColumnWise(" + name + ")", sum);
	}

	public int getLength()
	{
		return rowDimension * colDimension;
	}
	
	static private String significantFiguresString(double value,int sig)
	{
		MathContext mathContext = new MathContext(sig,RoundingMode.HALF_DOWN); 
		BigDecimal bigDecimal = new BigDecimal(value, mathContext); 
		String ret=bigDecimal.toEngineeringString();
		// RVK: Truncate the trailing zeroes to match UNIX output.
		if(ret.contains("."))
			while(ret.endsWith("0"))
				ret=ret.substring(0,ret.length()-1);
		return ret;
	}
	
	static public double significantFigures(double value, int sig)
	{
		if (Double.isNaN(value) || Double.isInfinite(value))
			return value;
		
		MathContext mathContext = new MathContext(sig,RoundingMode.HALF_DOWN); 
		BigDecimal bigDecimal = new BigDecimal(value, mathContext);
		return bigDecimal.doubleValue();
	}

	public void reducePrecision(int sig)
	{
		for (int i = 0; i < rowDimension; i++)
		{
			double[] row = data[i];
			for (int j = 0; j < colDimension; j++)
			{
				row[j]=significantFigures(row[j],sig);
			}
		}
	}
	
	public void writeAsText(String fileName) throws IOException
	{
		writeAsText(fileName, 18);
	}
	
	public void writeAsText(String fileName,int sig) throws IOException
	{
		System.out.println("Writing " + this + " to " + fileName);
		PrintStream ps = null;
		try
		{
			MultiProtocolFile mpf = new MultiProtocolFile(fileName);
			OutputStream os = new SafeMultiProtocolFileOutputStream(mpf);
			ps = new PrintStream(os);

			for (int i = 0; i < rowDimension; i++)
			{
				double[] row = data[i];
				StringBuffer sb = new StringBuffer();
				for (int j = 0; j < colDimension; j++)
				{
					if (j > 0)
						sb.append(" ");
					sb.append(significantFiguresString(row[j],sig));
					row[j]=significantFigures(row[j],sig);
				}
				sb.append('\n');
				ps.print(sb.toString());
			}
		}
		finally
		{
			if (ps != null)
				ps.close();
		}

		System.out.println("Done.");
	}

	public Thread writeArrayAsync(final String zipFilePath)
	{
		Thread worker = new Thread(new Runnable()
		{
			@Override
			public void run()
			{
				try
				{
					writeArray(zipFilePath);
				}
				catch (IOException e)
				{
					log.error("Failed to write " + zipFilePath, e);
				}
				finally
				{
					asyncWriteInProgress.decrementAndGet();
				}
			}
		});
		// No daemon! Process MUST wait on termination .....
		// worker.setDaemon(true);
		asyncWriteInProgress.incrementAndGet();
		worker.start();
		return worker;
	}
	

	public void writeAsCsvFile(String fileName, String delimiter) throws IOException
	{
		writeAsCsvFile(fileName, null, delimiter);
	}

	public void writeAsCsvFile(String fileName, List<String> header, String delimiter) throws IOException
	{
		MultiProtocolFile mpf = new MultiProtocolFile(fileName);
		OutputStream os = new SafeMultiProtocolFileOutputStream(mpf);
		PrintStream ps = new PrintStream(new BufferedOutputStream(os));

		if (header != null)
		{
			StringBuffer sb = new StringBuffer();
			for (String colName : header)
			{
				sb.append(colName);
				sb.append(delimiter);
			}
			ps.println(sb.toString());
		}

		for (int i = 0; i < data.length; i++)
		{
			StringBuffer sb = new StringBuffer();

			double[] row = data[i];
			for (int j = 0; j < row.length; j++)
			{
				if (j > 0)
					sb.append(delimiter);
				sb.append(String.valueOf(row[j]));
			}

			ps.println(sb.toString());

		}
		ps.close();

	}

	public void writeArray(String fileName) throws IOException
	{
//		long t0 = System.currentTimeMillis();
		System.out.println("Writing " + this + " to " + fileName);
		OutputStream os = null;
		try
		{
			MultiProtocolFile mpf = new MultiProtocolFile(fileName);
			os = new SafeMultiProtocolFileOutputStream(mpf);
			writeArray(os);
		}
		finally
		{
			if (os != null)
				os.close();
		}

//		long t1 = System.currentTimeMillis();
//		System.out.println("Done, " + (t1 - t0) + "ms");
	}
	
//	static public enum ComparisonMode
//	{
//		NONE,WRITE_TO_FILE,COMPARE_TO_FILE
//	};
//	static ComparisonMode comparisonMode=ComparisonMode.NONE;
//	static String comparePath="";
//	static public void setComparisonMode(ComparisonMode c)
//	{
//		comparisonMode=c;
//	}
//	static public void setComparisonPath(String p)
//	{
//		comparePath=p;
//	}
//	/**
//	 * Compare to the matrix stored in the stated file, in the current compare path (setComparisonPath)
//	 * but only if comparisonMode (setComparisonMode) is not NONE.
//	 */
//	public void compare(String fileName) 
//	{
//		if(comparisonMode==ComparisonMode.NONE)
//			return;
//		Mat compare;
//		try
//		{
//			compare = Mat.readArray("compare",comparePath+fileName,getRowDimension(),0);
//			if(!this.equals(compare))
//			{
//				System.err.print("Comparison failed for matrix "+name);
//			}
//		}
//		catch (InterruptedException e)
//		{
//			e.printStackTrace();
//		}
//		catch(IOException e)
//		{
//			System.err.print("Comparison failed for matrix "+name+" due to nonexistent file "+comparePath+fileName);
//			e.printStackTrace();
//		}
//	}
//	public void writeOrCompare(String fileName)
//	{
//		if(comparisonMode==ComparisonMode.WRITE_TO_FILE)
//		{
//	        File dir = new File(comparePath);
//	        dir.mkdirs();
//	        try
//	        {
//	        	writeArray(comparePath+fileName);
//	        }
//	        catch(IOException e)
//	        {
//				System.err.print("Failed to write matrix "+name+" to file "+comparePath+fileName);
//				e.printStackTrace();
//	        }
//		}
//		else if(comparisonMode==ComparisonMode.COMPARE_TO_FILE)
//		{
//			compare(fileName);
//		}
//			
//	}
		
		
	protected void writeDoubleArray(DataOutputStream dos, double[][] d, int j) throws IOException
	{
		ByteBuffer b = ByteBuffer.allocate(8 * d.length);
		
		if(outputFileFormat==FileFormat.DEFAULT)
		{
			// Copy the java array's contents to the native memory
			DoubleBuffer db = b.asDoubleBuffer();
	
			for (int i = 0; i < d.length; i++)
			{
				db.put(data[i][j]);
			}
		}
		else
		{
			// Copy the java array's contents to the native memory
			LongBuffer lb = b.asLongBuffer();
			
			for (int i = 0; i < d.length; i++)
			{
				lb.put( Long.reverseBytes( Double.doubleToLongBits(data[i][j]) ) );
			}
		}
		
		dos.write(b.array(), 0, b.array().length);
	}

	public void writeArray(OutputStream os) throws IOException
	{
		writeArray(os, 0, colDimension);
	}
	
	public void writeArray(OutputStream os, int start, int end) throws IOException
	{
		DataOutputStream dos = new DataOutputStream(new BufferedOutputStream(os));

		// RVK: replacing a lot of writeDouble's with a few writeDoubleArray's,
		// but the data is written in column-major order so the benefit is
		// limited by
		// the need to transpose.
		for (int j = start; j < end; j++)
		{
			writeDoubleArray(dos, data, j);
		}
		dos.flush();
	}

	public static Mat fillRange(String name, int start, int end)
	{
		double[] result = new double[end - start + 1];
		for (int i = 0; i < result.length; i++)
			result[i] = start + i;
		Mat m = new Mat(name, new double[][] { result });

		if (verbose)
			m.print();

		return m;
	}

	public static Mat mergeCols(String name, Mat... mv)
	{
		int rows = 0;
		int cols = 0;
		for (Mat m : mv)
		{
			if (m != null)
			{
				cols += m.colDimension;
				rows = m.rowDimension;
			}
		}

		double[][] result = new double[rows][cols];

		int ofs = 0;
		for (Mat m : mv)
		{
			if (m == null)
				continue;

			double[][] d = m.getData();
			int thisCols = m.colDimension;
			for (int j = 0; j < thisCols; j++)
			{
				for (int i = 0; i < rows; i++)
				{
					result[i][j + ofs] = d[i][j];
				}
			}
			ofs += m.colDimension;
		}

		return new Mat(name, result);
	}

	public static Mat mergeRows(String name, Mat... mv)
	{
		int cols = 0;
		int rows = 0;
		for (Mat m : mv)
		{
			if (m != null)
			{
				rows += m.rowDimension;
				cols = m.colDimension;
			}
		}

		double[][] result = new double[rows][cols];

		int ii = 0;
		for (Mat m : mv)
		{
			if (m == null)
				continue;

			double[][] d = m.getData();

			for (int i = 0; i < d.length; i++)
			{
				double[] a = d[i];
				double[] c = result[ii];
				for (int j = 0; j < cols; j++)
				{
					c[j] = a[j];
				}
				ii++;
			}
		}

		return new Mat(name, result);
	}
	
	// Fast merge. Stores back rows (NO COPY) from Mat m to this, starting with ofs
	public void setInnerRows(Mat m, int ofs)
	{
		double[][] mData = m.getData();
		for (int i = 0; i < mData.length; i++)
			data[i + ofs] = mData[i];
	}

	public static int[] asInteger(double[] array)
	{
		int[] result = new int[array.length];
		for (int i = 0; i < array.length; i++)
			result[i] = (int) array[i];
		return result;
	}

	public static int[] asInteger(List<Integer> list)
	{
		int[] result = new int[list.size()];
		for (int i = 0; i < result.length; i++)
			result[i] = list.get(i);
		return result;
	}

	public Mat(double[][] data)
	{
		this(concat("[", data.length, ",", data.length == 0 ? 0 : data[0].length, "]"), data);
	}

	public Mat(String name, double[][] data)
	{
		this.name = name;
		this.data = data;
		this.rowDimension = data.length;
		this.colDimension = data.length > 0 ? data[0].length : 0;
	}

	public Mat(String name, int rows, int cols)
	{
		this.name = name;
		this.data = new double[rows][cols];
		this.rowDimension = rows;
		this.colDimension = cols;
	}
	
		
	/**
	 * The next operation will be performed on this Mat itself....
	 * 
	 * @return
	 */
	public Mat inner()
	{
		isInner = true;
		return this;
	}
	
	
	/**
	 * Returns a new Mat in case isInner, otherwise returns this.
	 * isInner will be set to false in any case.
	 * 
	 * @param name
	 * @param array
	 * @return
	 */
	private Mat innerOrNew(String name, double[][] array)
	{		
		boolean wasInner = isInner;
		isInner = false;
		
		if (wasInner && data == array)
		{
			this.name = name;
			return this;
		}
		
		return new Mat(name, array);
	}
	
	private double[][] createResult()
	{
		return isInner ? data : new double[rowDimension][colDimension];
	}

	public String getName()
	{
		return name;
	}

	public void setName(String name)
	{
		this.name = name;
	}

	/**
	 * Fill rows from top to bottom.
	 * 
	 * @param values
	 */
	public void fillRows(double... values)
	{
		for (int i = 0; i < rowDimension; i++)
		{
			double[] row = data[i];
			int j = 0;
			for (double d : values)
				row[j++] = d;
		}
	}
	
	public void fill(double value)
	{
		for (int i = 0; i < rowDimension; i++)
		{
			double[] row = data[i];
			for (int j = 0; j < colDimension; j++)
				row[j] = value;
		}
	}

	/**
	 * Values is a column vector fill columns from
	 * left to right.
	 * 
	 * @param values
	 */
	public void fillColumns(double... values)
	{
		for (int j = 0; j < colDimension; j++)
		{
			int i = 0;
			for (double d : values)
				data[i++][j] = d;
		}
	}
	
	public static double[] asDoubleArray(float[] fa)
	{
		double[] da = new double[fa.length];
		for (int i = 0; i < fa.length; i++)
			da[i] = fa[i];
		return da;
	}

	
	public double[][] getData()
	{
		return data;
	}

	public Mat getColumnMatrix(int column)
	{
		return selectColumns(new int[] { column }, 0);
	}

	public Mat getRowMatrix(int row)
	{
		return selectRows(new int[] { row }, 0);
	}

	public Mat selectRows(Mat selector)
	{
		int[] rowArray;

		if (selector.rowDimension > 1)
			rowArray = asInteger(selector.getColumn(0));
		else
			rowArray = asInteger(selector.getRow(0));

		// !! in matrix objects we assume counting from 1:n !
		Mat m = selectRows(rowArray, -1);
		m.setName(concat(name, "(", selector.name, ",:)"));

		if (verbose)
			m.print();

		return m;
	}

	public Mat selectRows(int[] rowArray, int offset)
	{
		double[][] result = new double[rowArray.length][colDimension];

		for (int i = 0; i < rowArray.length; i++)
		{
			double[] row = result[i];
			for (int j = 0; j < colDimension; j++)
			{
				row[j] = data[rowArray[i] + offset][j];
			}
		}

		return new Mat(result);
	}
	
	public Mat dropRow(int rowIndex)
	{
		double[][] result = new double[rowDimension - 1][colDimension];

		int k = 0;
		for (int i = 0; i < rowDimension; i++)
		{
			if (i == rowIndex)
				continue;
			
			double[] row = result[k];
			for (int j = 0; j < colDimension; j++)
			{
				row[j] = data[i][j];
			}
			k++;
		}

		return new Mat(result);		
	}

	public Mat selectColumns(Mat selector)
	{
		int[] colArray;

		if (selector.rowDimension > 1)
			colArray = asInteger(selector.getColumn(0));
		else
			colArray = asInteger(selector.getRow(0));

		// !! in matrix objects we assume counting from 1:n !
		Mat m = selectColumns(colArray, -1);
		m.setName(concat(name, "(:,", selector.name, ")"));

		if (verbose)
			m.print();

		return m;
	}

	public Mat selectColumns(int[] colArray, int offset)
	{
		double[][] result = new double[rowDimension][colArray.length];

		for (int i = 0; i < rowDimension; i++)
		{
			double[] a = data[i];
			double[] c = result[i];
			for (int j = 0; j < colArray.length; j++)
			{
				c[j] = a[colArray[j] + offset];
			}
		}
		return new Mat(result);
	}

	public double get(int row, int col)
	{
		return data[row][col];
	}

	public void set(int row, int col, double value)
	{
		data[row][col] = value;
	}

	public void set(Mat mask, double v)
	{
		for (int i = 0; i < rowDimension; i++)
		{
			double[] a = data[i];
			double[] c = mask.data[i];
			for (int j = 0; j < colDimension; j++)
			{
				if (c[j] != 0.)
					a[j] = v;
			}
		}
	}

	public Mat setColumns(int[] colArray, int offset, Mat mat)
	{
		double[][] result = new double[rowDimension][colDimension];
		double[][] matData = mat.getData();

		// copy
		for (int i = 0; i < rowDimension; i++)
		{
			double[] a = data[i];
			double[] c = result[i];
			for (int j = 0; j < colDimension; j++)
			{
				c[j] = a[j];
			}
		}

		for (int k = 0; k < colArray.length; k++)
		{
			int idx = colArray[k] + offset;

			for (int i = 0; i < rowDimension; i++)
			{
				result[i][idx] = matData[i][k];
			}
		}

		return new Mat(result);
	}

	public void setInnerColumns(int[] colArray, int offset, Mat mat)
	{
		double[][] matData = mat.getData();
		for (int k = 0; k < colArray.length; k++)
		{
			int idx = colArray[k] + offset;
			for (int i = 0; i < rowDimension; i++)
				data[i][idx] = matData[i][k];
		}
	}

	public double[] getRow(int i)
	{
		double[] a = data[i];
		int n = a.length;

		double[] d = new double[n];
		for (int j = 0; j < n; j++)
			d[j] = a[j];

		return d;
	}

	public double[] getColumn(int col)
	{
		double[] d = new double[data.length];
		for (int i = 0; i < data.length; i++)
			d[i] = data[i][col];

		return d;
	}
		

	public Mat dotDivide(Mat m)
	{
		double[][] d = m.getData();

		int rows = m.rowDimension;
		int cols = m.colDimension;
		double[][] result = createResult();

		for (int i = 0; i < rows; i++)
		{
			double[] a = data[i];
			double[] b = d[i];
			double[] c = result[i];
			for (int j = 0; j < cols; j++)
				c[j] = a[j] / b[j];
		}

		Mat res = innerOrNew(concat("(", name, "./", m.name, ")"), result);

		if (verbose)
			res.print();

		return res;
	}

	public Mat dotMultiply(Mat m)
	{
		final double[][] result = createResult();

		for (int i = 0; i < rowDimension; i++)
		{
			double[] a = data[i];
			double[] b = m.data[i];
			double[] c = result[i];
			for (int j = 0; j < colDimension; j++)
				c[j] = b[j] * a[j];
		}

		return innerOrNew(concat("(", this.name, ".*", m.name, ")"), result);
	}

	public Mat dotPow(double exponent)
	{
		final double[][] result = createResult();

		for (int i = 0; i < rowDimension; i++)
		{
			double[] a = data[i];
			double[] c = result[i];
			for (int j = 0; j < colDimension; j++)
				c[j] = Math.pow(a[j], exponent);
		}

		return innerOrNew(concat("(", this.name, ".^", exponent, ")"), result);
	}

	public Mat dotReciprocal(double multiplier)
	{
		final double[][] result = createResult();

		for (int i = 0; i < rowDimension; i++)
		{
			double[] a = data[i];
			double[] c = result[i];
			for (int j = 0; j < colDimension; j++)
				c[j] = multiplier / a[j];
		}

		return innerOrNew(concat("(", multiplier, "./", name, ")"), result);
	}

	public Mat sqrt()
	{
		final double[][] result = createResult();
		
		MatRowWorker job = new MatRowWorker(concat("sqrt(", name, ")"), result, data, null) {

			public int calculationsPerThread()
			{
				return 10000;
			}

			public void run()
			{
				for (int i = iStart; i < iEnd; i++)
				{
					double[] a = A[i];
					double[] c = result[i];
					for (int j = 0; j < a.length; j++)
					{
						double v = Math.sqrt(a[j]);
						if (Double.isNaN(v))
						{
							if (verbose)
								System.out.println("[" + i + "," + j + "] = NaN, setting to 0.0");
							v = 0;
						}
						c[j] = v;
					}
				}			
			}
		};
		return executeMatRowWorkers(job, simpleMethodMaxThreads);
	}

	public Mat exp()
	{
		final double[][] result = createResult();

		MatRowWorker job = new MatRowWorker(concat("exp(", name, ")"), result, data, null) {

			public void run()
			{
				for (int i = iStart; i < iEnd; i++)
				{
					double[] a = A[i];
					double[] c = result[i];
					for (int j = 0; j < a.length; j++)
					{
						double val = Math.exp(a[j]);
						if (Double.isNaN(val))
						{
							if (verbose)
								System.out.println("[" + i + "," + j + "] = NaN, setting to 0.0");
							val = 0;
						}
						c[j] = val;
					}
				}			
			}
		};
		return executeMatRowWorkers(job, simpleMethodMaxThreads);
	}
	
	public Mat ln()
	{
		final double[][] result = createResult();

		MatRowWorker job = new MatRowWorker(concat("exp(", name, ")"), result, data, null) {

			public void run()
			{
				for (int i = iStart; i < iEnd; i++)
				{
					double[] a = A[i];
					double[] c = result[i];
					for (int j = 0; j < a.length; j++)
					{
						c[j] = Math.log(a[j]);
					}
				}			
			}
		};
		return executeMatRowWorkers(job, simpleMethodMaxThreads);
	}

	public Mat vectorMultiply(Mat v)
	{
		// v must have one col, rowDimensions equal!
		
		double[][] result = createResult();
		for (int i = 0; i < rowDimension; i++)
		{
			double[] a = data[i];
			double[] c = result[i];
			double f = v.data[i][0];
			for (int j = 0; j < colDimension; j++)
			{
				c[j] = a[j] * f;
			}
		}
		return innerOrNew(concat("(", name, "*", v, ")"), result);
	}

	public Mat scalarMultiply(double f)
	{
		double[][] result = createResult();
		for (int i = 0; i < rowDimension; i++)
		{
			double[] a = data[i];
			double[] c = result[i];
			for (int j = 0; j < colDimension; j++)
			{
				c[j] = a[j] * f;
			}
		}
		return innerOrNew(concat("(", name, "*", f, ")"), result);
	}

	public Mat scalarAdd(double d)
	{
		double[][] result = createResult();

		for (int i = 0; i < rowDimension; i++)
		{
			double[] a = data[i];
			double[] c = result[i];
			for (int j = 0; j < colDimension; j++)
			{
				c[j] = a[j] + d;
			}
		}

		return innerOrNew(concat("(", name, "+", d, ")"), result);
	}
	
	
	public Mat scalarSubstract(double d)
	{
		double[][] result = createResult();

		for (int i = 0; i < rowDimension; i++)
		{
			double[] a = data[i];
			double[] c = result[i];
			for (int j = 0; j < colDimension; j++)
			{
				c[j] = a[j] - d;
			}
		}

		return innerOrNew(concat("(", name, "-", d, ")"), result);
	}


	// TODO
	public Mat elimnans(double substitute4NaN)
	{
		for (int i = 0; i < rowDimension; i++)
		{
			double[] a = data[i];
			for (int j = 0; j < colDimension; j++)
			{
				if (Double.isNaN(a[j]))
				{
					System.err.println("[" + i + "," + j + "] = NaN, set to " + substitute4NaN);
					a[j] = substitute4NaN;
				}
			}
		}
		return this;
	}

	public Mat add(Mat m)
	{
		double[][] result = createResult();
		double[][] mData = m.getData();

		for (int i = 0; i < rowDimension; i++)
		{
			double[] a = data[i];
			double[] b = mData[i];
			double[] c = result[i];
			for (int j = 0; j < colDimension; j++)
			{
				c[j] = a[j] + b[j];
			}
		}
		return innerOrNew(concat("(", name, "+", m.name, ")"), result);
	}

	public Mat substract(Mat m)
	{
		double[][] result = createResult();
		double[][] mData = m.getData();

		for (int i = 0; i < rowDimension; i++)
		{
			double[] a = data[i];
			double[] b = mData[i];
			double[] c = result[i];
			for (int j = 0; j < colDimension; j++)
			{
				c[j] = a[j] - b[j];
			}
		}
		return innerOrNew(concat("(", name, "-", m.name, ")"), result);
	}
	
	
	private static void postMultiplyByDiagonals(double[][] A, double[] diag, int iStart, int iEnd)
	{
		final int cols = diag.length;
		for (int i = iStart; i < iEnd; i++)
		{
			double[] rowA = A[i];
			for (int j = 0; j < cols; j++)
				rowA[j] *= diag[j];
		}
	}

	public void postMultiplyByDiagonals(double [] diag)
	{
		postMultiplyByDiagonals(data,diag,0,data.length);
	}
	
	
	private static void preMultiplyByDiagonals(double[][] A, double[] diag, int iStart, int iEnd)
	{
		final int cols = A[0].length;
		for (int i = iStart; i < iEnd; i++)
		{
			double[] rowA = A[i];
			for (int j = 0; j < cols; j++)
				rowA[j] *= diag[i];
		}
	}

	public void preMultiplyByDiagonals(double [] diag)
	{
		preMultiplyByDiagonals(data,diag,0,data.length);
	}
	
	
	public Mat multiplyByTranspose(Mat m)
	{
		double[][] A = data;
		double[][] B = m.data;
		MatRowWorker job = new MatRowWorker(concat("(", name, "*", m.name, "')"), new double[A.length][B.length], A, B) {
			public void run() 
			{
				multiplyByTranspose(result, A, B, iStart, iEnd);				
			}
		};
		return executeMatRowWorkers(job);
	}
	
	
	public Mat multiply(Mat m)
	{
		double[][] A = data;
		double[][] B = m.data;
		final int cols = B[0].length;
		
		MatRowWorker job = new MatRowWorker(concat("(", name, "*", m.name, ")"),new double[A.length][cols], A, B) {
			public void run()
			{
				multiply(result, A, B, iStart, iEnd);				
			}
		};
		Mat resultMat =  executeMatRowWorkers(job);


		// TEST
		if (testScaling)
		{
			System.out.println("verifying " + name + "*" + m.name + "...");
			final double[][] resTest = new double[A.length][cols];
			multiply(resTest, A, B, 0, A.length);
			if (!isEqual(resTest, resultMat.data, 0.))
				throw new RuntimeException("MultipleWorker multiply() failed");
		}

		return resultMat;
	}
	
	
	public boolean isSymmetric()
	{
		if (rowDimension != colDimension)
			return false;
		
		int n = rowDimension;
		
		for (int i = 0; i < n; i++)
		{
			double[] rowi = data[i];
			for (int j = 0; j < n; j++)
			{
				if (rowi[j] != data[j][i])
					return false;
			}
		}
				
		return true;
	}
	
	public boolean isDiagonal()
	{
		return isDiagonal(0.);
	}
	
	public boolean isDiagonal(double threshold)
	{
		for (int i = 0; i < rowDimension; i++)
		{
			double[] rowi = data[i];
			for (int j = 0; j < colDimension; j++)
			{
				if (i == j)
					continue;
				
				double d = rowi[j];
				if (Math.abs(d) < threshold)
					d = 0.;
				
				if (d != 0.)
					return false;
			}
		}
				
		return true;
	}
	
	
	/**
	 * Check if any values is NaN
	 * @return
	 */
	public boolean isNaN()
	{
		for (int i = 0; i < rowDimension; i++)
		{
			double[] row = data[i];
			for (int j = 0; j < colDimension; j++)
			{
				 if(new Double(row[j]).isNaN())
					 return true;
			}
		}
		return false;
	}

	private Mat executeMatRowWorkers(MatRowWorker job)
	{
		return executeMatRowWorkers(job, getMaxThreads());
	}
	
	private Mat executeMatRowWorkers(MatRowWorker job, int proposedNumberOfThreads)
	{
		int maxThreads = Math.max(getMaxThreads(), proposedNumberOfThreads);
		
		double[][] A = job.getA();
		int numThreads = Math.min(maxThreads, job.proposeNumberOfThreads());
		int deltaRowsA = numThreads == 0 ? 0 : A.length / numThreads;
		if (executorService == null || numThreads < 2 || deltaRowsA < 2)
		{
			// single thread model
			job.setBounds(0, A.length);
			job.setIndex(0);
			job.run();
		}
		else
		{
			int iStart = 0;
			int iEnd = 0;
			List<Future<?>> futures = new LinkedList<Future<?>>();
			for (int t = 0; t < numThreads; t++)
			{
				iEnd = (t == numThreads - 1) ? A.length : iStart + deltaRowsA;
				MatRowWorker clone = job.shallowCopy();
				clone.setBounds(iStart, iEnd);
				clone.setIndex(t);
				futures.add(executorService.submit(clone));
				iStart = iEnd;
			}

			ExecutorUtils.getAllAndRemove(futures);			
		}
		
		double[][] r = job.getResult();
		return r == null ? null : innerOrNew(job.getName(), r);		
	}
	
	
	public boolean isEqual(Mat that)
	{
		return isEqual(that, 0.005);
	}

	public boolean isEqual(Mat that, double deviation)
	{
		return (isEqual(this.data, that.data, deviation));
	}

	private static boolean isEqual(double[][] a, double[][] b, double deviation)
	{
		if (a.length != b.length || a[0].length != b[0].length)
			return false;

		int rows = a.length;
		int cols = a[0].length;

		for (int i = 0; i < rows; i++)
		{
			double[] aRow = a[i];
			double[] bRow = b[i];
			for (int j = 0; j < cols; j++)
			{
				double v0 = aRow[j];
				double v1 = bRow[j];
				if (deviation != 0.)
				{
					double rel = Math.abs((v1 - v0) / (v0 == 0 ? 1 : v0));
					if (rel > deviation)
					{
						System.out.println(v0 + " != " + v1 + " at " + i + "/" + j);
						return false;
					}
				}
				else
				{
					if (v0 != v1)
						return false;
				}
			}
		}

		return true;
	}
/* RVK: This experimental OpenCL support is disabled.
	static CLContext context = null;
	static CLQueue queue = null;
	static CLKernel matrixMulKernel = null;
	static CLProgram program = null;
	private static synchronized void createCLContextAndProgram()
	{
		// CLPlatform.DeviceFeature[]
		// preferredFeatures={CLPlatform.DeviceFeature.DoubleSupport};
		context = JavaCL.createBestContext();
		queue = context.createDefaultQueue();
		System.out.println("multiplyCL created CLContext...");
		// Read the program sources and compile them :
		String src = null;
		try
		{
			src = IOUtils.readText(JavaCLTest.class.getResource("Mat.cl"));
		}
		catch (IOException e)
		{
			e.printStackTrace();
		}
		program = context.createProgram(src);
		matrixMulKernel = program.createKernel("matrixMul");
		System.out.println("multiplyCL created CLProgram...");
	}
	public static synchronized void multiplyCL(double[][] res, double[][] A, double[][] B)
	{
		System.out.println("multiplyCL " + A.length + "x" + A[0].length + "*" + B.length + "x"
				+ B[0].length + "...");
		if (context == null)
		{
			createCLContextAndProgram();
		}
		final int ROW = 1024;
		final int COL = 1024;
		final int wA = A[0].length;
		final int hB = B.length;
		int memSizeB = hB * COL;

		Pointer<Float> bPtr = allocateFloats(memSizeB);
		// Create OpenCL input buffers (using the native memory pointers aPtr
		// and bPtr) :

		int memSizeC = ROW * COL;
		// Create an OpenCL output buffer:
		CLBuffer<Float> out = context.createFloatBuffer(Usage.Output, memSizeC);
		for (int row = 0; row < A.length; row += ROW)
		{
			int hA = java.lang.Math.min(ROW, A.length - row);
			int memSizeA = hA * wA;
			Pointer<Float> aPtr = allocateFloats(memSizeA);

			for (int i = 0; i < hA; i++)
			{
				for (int j = 0; j < wA; j++)
					aPtr.set(i * wA + j, (float) A[row + i][j]);
			}
			CLBuffer<Float> a = context.createFloatBuffer(Usage.Input, aPtr, true);

			for (int col = 0; col < B[0].length; col += COL)
			{
				int wB = java.lang.Math.min(COL, B[0].length - col);
				for (int i = 0; i < hB; i++)
				{
					for (int j = 0; j < wB; j++)
						bPtr.set(j * hB + i, (float) B[i][col + j]);
					// We've converted B to column-major for speed here.
				}
				CLBuffer<Float> b = context.createFloatBuffer(Usage.Input, bPtr, true);
				matrixMulKernel.setArgs(a, b, out, wA, wB);
				// System.out.append("(a,b,c,wA="+wA+",wB="+wB+"){"+hA+","+wB+"}");
				CLEvent addEvt = matrixMulKernel.enqueueNDRange(queue, new int[] { hA, wB });
				// System.out.append('.');
				Pointer<Float> outPtr = out.read(queue, addEvt);
				// blocks until finished
				for (int i = 0; i < hA; i++)
				{
					for (int j = 0; j < wB; j++)
						res[row + i][col + j] = outPtr.get(i * wB + j);
				}
			}
			// System.out.println("");
		}
		// System.out.println("Done");
	}*/

	private static void multiply(double[][] res, double[][] A, double[][] B, int iStart, int iEnd)
	{
		final int rowsB = B.length;
		final int cols = B[0].length;

		for (int i = iStart; i < iEnd; i++)
		{
			double[] row = res[i];
			double[] rowA = A[i];
			for (int j = 0; j < cols; j++)
			{
				double c = row[j];			
				for (int k = 0; k < rowsB; k++)
					c += rowA[k] * B[k][j];
				row[j] = c;
			}
		}
	}
	
	private static void multiplyByTranspose(double[][] res, double[][] A, double[][] B, int iStart, int iEnd)
	{
		final int rowsB = B[0].length;
		final int cols = B.length;

		for (int i = iStart; i < iEnd; i++)
		{
			double[] row = res[i];
			double[] rowA = A[i];
			for (int j = 0; j < cols; j++)
			{
				double[] colB = B[j];
				double c = 0.0;
				for (int k = 0; k < rowsB; k++)
					c += rowA[k] * colB[k];
				row[j] = c;
			}
		}
	}
	
	public Mat transpose()
	{
		final double[][] result = new double[colDimension][rowDimension];

		for (int i = 0; i < rowDimension; i++)
		{
			double[] a = data[i];
			for (int j = 0; j < colDimension; j++)
			{
				result[j][i] = a[j];
			}
		}

		return new Mat(name + "'", result);
	}

	public Mat log()
	{
		double[][] r = createResult();
		
		MatRowWorker job = new MatRowWorker(concat("log(", name, ")"), r, data, null) {

			public int calculationsPerThread()
			{
				return 10000;
			}

			public void run()
			{
				for (int i = iStart; i < iEnd; i++)
				{
					double[] a = A[i];
					double[] c = result[i];
					for (int j = 0; j < a.length; j++)
					{
						c[j] = Math.log(a[j]);
					}
				}			
			}
		};
		return executeMatRowWorkers(job, simpleMethodMaxThreads);

		
	}

	/**
	 * Variance by column i.e. colums are seen as row vectors
	 * @return
	 */
	public Mat var()
	{
		double w = 1. / (rowDimension - 1);
		double wMean = 1. / (rowDimension);	// Bugfix, formerly rowDimension - 1

		double[] mean = new double[colDimension];
		double[] var = new double[colDimension];

		for (int i = 0; i < rowDimension; i++)
		{
			double[] row = data[i];
			for (int j = 0; j < colDimension; j++)
				mean[j] += row[j] * wMean;
		}

		for (int i = 0; i < rowDimension; i++)
		{
			double[] row = data[i];
			for (int j = 0; j < colDimension; j++)
			{
				double d = row[j] - mean[j];
				var[j] += d * d * w;
			}
		}

		return new Mat(concat("var(", name, ")"), new double[][] { var });
	}
	
	/**
	 * Covariance of columns
	 * @return
	 */
//	public Mat cov()
//	{
//		// Ai - mean
//		Mat A = this.copy(name);
//		Mat mean = mean();
//		
//		double[][] d = A.data;
//		double[] md = mean.data[0];
//		for (int i = 0; i < A.rowDimension; i++)
//		{
//			double[] rows = d[i];
//			for (int j = 0; j < A.colDimension; j++)
//				rows[j] -= md[j];
//		}
//		
//		// E[A' x A]
//		Mat cov = A.transpose().multiply(A);
//		cov.inner().scalarMultiply(1. / (A.rowDimension - 1));
//		cov.setName(concat("cov(", name, ")"));
//		return cov;
//	}
	
	public Mat cov()
	{
		Mat A = distance2Mean(this);
		
		// E[A' x A]
		Mat cov = A.transpose().multiply(A);
		cov.inner().scalarMultiply(1. / (A.rowDimension - 1));
		cov.setName(concat("cov(", name, ")"));
		return cov;
	}

	public Mat cov(Mat B)
	{
		Mat tmpA = distance2Mean(this);
		Mat tmpB = distance2Mean(B);
		
		Mat cov = tmpA.transpose().multiply(tmpB);
		cov.inner().scalarMultiply(1. / (rowDimension - 1));
		cov.setName(concat("cov(", name, ",", B.name, ")"));
		
		return cov;		
	}

	
//	public static Mat cov(Mat A, Mat B)
//	{
//		Mat tmpA = distance2Mean(A);
//		Mat tmpB = distance2Mean(B);
//		
//		Mat cov = tmpA.transpose().multiply(tmpB);
//		cov.inner().scalarMultiply(1. / (A.rowDimension - 1));
//		
//		return cov;		
//	}
	
	private static Mat distance2Mean(Mat m)
	{
		Mat dist = m.copy("dist");
		Mat mean = m.mean();

		double[][] d = dist.data;
		double[] md = mean.data[0];
		for (int i = 0; i < m.rowDimension; i++)
		{
			double[] rows = d[i];
			for (int j = 0; j < m.colDimension; j++)
				rows[j] -= md[j];
		}
		
		return dist;
	}
	
	public Mat correlation(Mat B)
	{
		Mat cor = cov(B);
		Mat varA = var();
		Mat varB = B.var();
		double[] varARow = varA.data[0];
		double[] varBRow = varB.data[0];
	
		int n = cor.rowDimension;
		int m = cor.colDimension;
		for (int i = 0; i < n; i++)
		{
			double[] row = cor.data[i];
			for (int j = 0; j < m; j++)
			{
				double a = row[j];
				double b = varARow[i] * varBRow[j];
				
				if (b == 0.)
					row[j] = 0.;
				else
					row[j] = a / Math.sqrt( b );
			}
		}
		
		cor.setName(concat("correlation(", name, ")"));
		return cor;
	}
	
	
	public Mat correlation()
	{
		Mat cor = cov();
		Mat var = var();
		double[] varRow = var.data[0];
		
		int n = cor.rowDimension;
		for (int i = 0; i < n; i++)
		{
			double[] row = cor.data[i];
			for (int j = 0; j < n; j++)
			{
				double a = row[j];
				double b = varRow[i] * varRow[j];
				
				if (a == 0. && b == 0.)
					row[j] = i == j ? 1 : 0;
				else
					row[j] = a / Math.sqrt( b );
			}
		}
		
		cor.setName(concat("correlation(", name, ")"));
		return cor;
	}

	public void print(String name)
	{
		print(name, 5, 10);
	}

	public void print()
	{
		print(this.name, 5, 10);
	}
	
	public void dump(PrintStream ps, String name)
	{
		if (ps != null)
			print(ps, name, getRowDimension(), getColumnDimension());
	}

	public void print(String name, int rows, int cols)
	{
		print(System.out, name, rows, cols);
	}
	
	public void print(PrintStream ps, String name, int rows, int cols)
	{
		rows = Math.min(rows, data.length);
		cols = Math.min(cols, data[0].length);

		ps.println((name == null ? this.name : name) + ": [" + rowDimension + ", " + colDimension + "]");

		for (int i = 0; i < rows; i++)
		{
			StringBuffer sb = new StringBuffer();
			for (int j = 0; j < cols; j++)
			{
				if (j > 0)
					sb.append(", ");
				sb.append(new Formatter().format("%.7f", data[i][j]));
			}
			ps.println(sb);
		}
	}

	public int getRowDimension()
	{
//		if (data.length != this.rowDimension)
//			throw new RuntimeException("This is very bad ....");
		
//		return data.length;
		return rowDimension;
	}

	public int getColumnDimension()
	{
//		if ((data.length > 0 ? data[0].length : 0) != this.colDimension)
//			throw new RuntimeException("This is very bad ....");
//
//		return data.length > 0 ? data[0].length : 0;
		return colDimension;
	}

	public String toString()
	{
		return concat(name, " [", rowDimension, ", ", colDimension, "]");

	}

	public Mat copy(String name)
	{
		double[][] clone = new double[rowDimension][colDimension];
		for (int i = 0; i < data.length; i++)
		{
			double[] src = data[i];
			double[] tgt = clone[i];
			for (int j = 0; j < src.length; j++)
				tgt[j] = src[j];
		}

		return new Mat(name, clone);
	}

	// linint2 can interpolate more than one value at once!
	
	// x = 1 by n list of x values
	// y = m by n list of y values
	// xnew = m by h list of new x values. A new y is to be found for each x, resulting in
	// ynew = m by h list of new y values.
	public static Mat linint2(String name, Mat x, Mat y, Mat xnew)
	{

		int len = Math.max(x.rowDimension, x.colDimension);
		if (len == 1)
		{
			// ynew=y(:,1)*ones(1,length(xnew));
			return y.selectColumns(new int[] { 0 }, 0).multiply(ones(1, len));
		}

		double[] xRow = x.getRow(0);
		double[][] yData = y.getData();
		double xMax = Double.MIN_VALUE;
		double xMin = Double.MAX_VALUE;
		int xMaxIdx = 0;
		int xMinIdx = 0;
		for (int i = 0; i < xRow.length; i++)
		{
			double v = xRow[i];
			if (v > xMax)
			{
				xMax = v;
				xMaxIdx = i;
			}
			if (v < xMin)
			{
				xMin = v;
				xMinIdx = i;
			}
		}

		double tmi;
		double tpl;
		double ymi[] = new double[yData.length];
		double ypl[] = new double[yData.length];

		// for t = xnew
		double[] xnewRow = xnew.getRow(0);

		// ynew is m by h
		double ynew[][] = new double[yData.length][xnewRow.length];

		for (int i = 0; i < xnewRow.length; i++)
		{
			double t = xnewRow[i];
			if (t >= xMax)
			{
				// tmi=min(x);
				// tpl=max(x);
				// ymi=y(:,find(x==tpl));
				// ypl=y(:,find(x==tpl));

				tmi = xMin;
				tpl = xMax;
				for(int j=0;j<yData.length;j++)
				{
					ymi [j]= yData[j][xMaxIdx];
					ypl [j]= yData[j][xMaxIdx];
				}
			}
			else if (t <= xMin)
			{
				// tmi=min(x);
				// tpl=max(x);
				// ymi=y(:,find(x==tmi));
				// ypl=y(:,find(x==tmi));

				tmi = xMin;
				tpl = xMax;
				for(int j=0;j<yData.length;j++)
				{
					ymi [j]= yData[j][xMinIdx];
					ypl [j]= yData[j][xMinIdx];
				}
			}
			else
			{
				// tmi=max((x<t).*x + (x>=t)*min(x));
				// tpl=min((x>=t).*x + (x<t)*max(x));
				// ypl=y(:,find(x==tpl));
				tmi = (x.lt(t).dotMultiply(x).add(x.gteq(t).inner().scalarMultiply(xMin))).max();
				tpl = (x.gteq(t).dotMultiply(x).add(x.lt(t).inner().scalarMultiply(xMax))).min();

				// ymi=y(:,find(x==tmi));
				for (int j = 0; j < xRow.length; j++)
				{
					double d = xRow[j];
					if (d == tmi)
					{
						for(int k=0;k<yData.length;k++)
							ymi[k]= yData[k][j];
						break;
					}
				}

				// ypl=y(:,find(x==tpl));
				for (int j = 0; j < xRow.length; j++)
				{
					double d = xRow[j];
					if (d == tpl)
					{
						for(int k=0;k<yData.length;k++)
							ypl[k]= yData[k][j];
						break;
					}
				}

			}

			// ynew=[ynew ((tpl-t)*ymi + (t-tmi)*ypl)/(tpl-tmi)];
			if(tpl == tmi)
			{
				for(int k=0;k<yData.length;k++)
					if(ymi[k]!=ypl[k])
						System.err.print("Can't interpolate "+name+" row "+i+" because 2 x's are the same but the y's are not!");
			}
			for(int k=0;k<yData.length;k++)
				ynew[k][i] = ((tpl - t) * ymi [k]+ (t - tmi) * ypl[k]) / (tpl - tmi);
		}

		return new Mat("ynew", ynew );
	}

	public static Mat linint(String name, Mat x, Mat y, Mat xnew)
	{

		int len = Math.max(x.rowDimension, x.colDimension);
		if (len == 1)
		{
			// ynew=y(:,1)*ones(1,length(xnew));
			return y.selectColumns(new int[] { 0 }, 0).multiply(ones(1, len));
		}

		double[] xRow = x.getRow(0);
		double[] yRow = y.getRow(0);
		double xMax = Double.MIN_VALUE;
		double xMin = Double.MAX_VALUE;
		int xMaxIdx = 0;
		int xMinIdx = 0;
		for (int i = 0; i < xRow.length; i++)
		{
			double v = xRow[i];
			if (v > xMax)
			{
				xMax = v;
				xMaxIdx = i;
			}
			if (v < xMin)
			{
				xMin = v;
				xMinIdx = i;
			}
		}

		double tmi;
		double tpl;
		double ymi = 0;
		double ypl = 0;;

		// for t = xnew
		double[] xnewRow = xnew.getRow(0);

		double ynew[] = new double[xnewRow.length];

		for (int i = 0; i < xnewRow.length; i++)
		{
			double t = xnewRow[i];
			if (t >= xMax)
			{
				// tmi=min(x);
				// tpl=max(x);
				// ymi=y(:,find(x==tpl));
				// ypl=y(:,find(x==tpl));

				tmi = xMin;
				tpl = xMax;
				ymi = yRow[xMaxIdx];
				ypl = yRow[xMaxIdx];
			}
			else if (t <= xMin)
			{
				// tmi=min(x);
				// tpl=max(x);
				// ymi=y(:,find(x==tmi));
				// ypl=y(:,find(x==tmi));

				tmi = xMin;
				tpl = xMax;
				ymi = yRow[xMinIdx];
				ypl = yRow[xMinIdx];
			}
			else
			{
				// tmi=max((x<t).*x + (x>=t)*min(x));
				// tpl=min((x>=t).*x + (x<t)*max(x));
				// ypl=y(:,find(x==tpl));
				tmi = (x.lt(t).dotMultiply(x).add(x.gteq(t).inner().scalarMultiply(xMin))).max();
				tpl = (x.gteq(t).dotMultiply(x).add(x.lt(t).inner().scalarMultiply(xMax))).min();

				// ymi=y(:,find(x==tmi));
				for (int j = 0; j < xRow.length; j++)
				{
					double d = xRow[j];
					if (d == tmi)
					{
						ymi = yRow[j];
						break;
					}
				}

				// ypl=y(:,find(x==tpl));
				for (int j = 0; j < xRow.length; j++)
				{
					double d = xRow[j];
					if (d == tpl)
					{
						ypl = yRow[j];
						break;
					}
				}

			}

			// ynew=[ynew ((tpl-t)*ymi + (t-tmi)*ypl)/(tpl-tmi)];
			ynew[i] = ((tpl - t) * ymi + (t - tmi) * ypl) / (tpl - tmi);
		}

		return new Mat("ynew", new double[][] { ynew });
	}
	
	/**
	 * Returns a Mat with one row containing the means of "this" column values.
	 * @return
	 */
	public Mat mean()
	{
		double[] mean = new double[colDimension];
		Mat meanMat = new Mat(concat("mean(", name, ")"), new double[][] { mean });
		
		double w = 1. / (rowDimension);

		for (int i = 0; i < rowDimension; i++)
		{
			double[] row = data[i];
			for (int j = 0; j < colDimension; j++)
				mean[j] += row[j] * w;
		}

		return meanMat;
	}

	/**
	 * return arithmetic mean of all values in matrix
	 * 
	 * @param x
	 * @return
	 */
	public double meanAllMatrix()
	{
		double mean = 0;

		for (int i = 0; i < rowDimension; i++)
		{
			double[] row = data[i];
			for (int j = 0; j < row.length; j++)
			{
				mean += row[j];
			}
		}
		mean /= (rowDimension * colDimension);

		return mean;
	}

	/**
	 * return standard deviation of all values in matrix
	 * 
	 * @return
	 */
	public double standardDevAllMatrix()
	{
		return standardDevAllMatrix(meanAllMatrix());
	}

	/**
	 * return standard deviation of all values in matrix
	 * 
	 * @param mean
	 * @return
	 */
	public double standardDevAllMatrix(double mean)
	{
		double standardDev = 0;
		for (int i = 0; i < rowDimension; i++)
		{
			double[] row = data[i];
			for (int j = 0; j < row.length; j++)
			{
				double d = mean - row[j];
				standardDev += (d*d);
			}
		}
		standardDev = Math.sqrt((standardDev) / (rowDimension * colDimension - 1));

		return standardDev;
	}

	/**
	 * return skewness of all values in matrix
	 * 
	 * @param d
	 * @return
	 */
	public double skewnessAllMatrix(double mean, double stdv)
	{
		double n = (rowDimension * colDimension);
		double sum3 = 0;
		for (int i = 0; i < rowDimension; i++)
		{
			double[] row = data[i];
			for (int j = 0; j < row.length; j++)
			{
				double d = row[j] - mean;
				sum3 += d*d*d;
			}
		}

		double skewness =sum3/(n*stdv*stdv*stdv);

		return skewness;
	}
	

	/**
	 * return skewness of all values in matrix
	 * 
	 * @return
	 */
	public double skewnessAllMatrix()
	{
		return skewnessAllMatrix(meanAllMatrix(),standardDevAllMatrix());
	}
	
	public double kurtosisAllMatrix()
	{
		return kurtosisAllMatrix(meanAllMatrix(),standardDevAllMatrix());
	}

	/**
	 * return kurtosis of all values in matrix
	 * 
	 * @param d
	 * @return
	 */
	public double kurtosisAllMatrix(final double mean, final double stdv)
	{
		double n = (rowDimension * colDimension);
		double sum4 = 0;
		for (int i = 0; i < rowDimension; i++)
		{
			double[] row = data[i];
			for (int j = 0; j < row.length; j++)
			{
				double d = row[j] - mean;
				sum4 +=d*d*d*d;
			}
		}

		double kurtosis = sum4/(n*stdv*stdv*stdv*stdv)-3;

		return kurtosis;
	}
	
	/**
	 * Returns a Mat with one row containing the kurtosis of "this" column values. 
	 * @return
	 */
	public Mat kurtosis(){
		double[] kurtosis = new double[colDimension];
		Mat d = this.transpose();
		for(int i=0;i<colDimension;i++){
			kurtosis[i]=d.selectRows(new int[]{i}, 0).kurtosisAllMatrix();
		}
		return new Mat(new double[][]{kurtosis});
	}

	/**
	 * Returns a Mat with one row containing the skewness of "this" column values. 
	 * @return
	 */
	public Mat skewness(){

		double[] skewness = new double[colDimension];
		Mat d = this.transpose();
		for(int i=0;i<colDimension;i++){
			skewness[i]=d.selectRows(new int[]{i}, 0).skewnessAllMatrix();
		}
		return new Mat(new double[][]{skewness});
	}
	

	/**
	 * Returns a Mat with one row containing the standardDev of "this" column values. 
	 * @return
	 */
	public Mat standardDev(){
		double[] standardDev = new double[colDimension];
		Mat d = this.transpose();
		for(int i=0;i<colDimension;i++){
			standardDev[i]=d.selectRows(new int[]{i}, 0).standardDevAllMatrix();
		}
		return new Mat(new double[][]{standardDev});
	}
	

	/**
	 * Generate a matrix using QR Decomposition using JAMA libraries
	 * 
	 * @param XX
	 * @return
	 */
	public Mat qr()
	{
		Mat QR = null;
		Matrix xxM = new Matrix(getData());
		QRDecomposition qrDec = xxM.qr();
		double[][] RR = qrDec.getR().getArray();
		double[][] HH = qrDec.getH().getArray();
		QR = triangularJoin(RR, HH);
		QR.setName("QR(" + name + ")");
		return QR;
	}

	/**
	 * Return the effective numerical rank, obtained from SVD
	 * 
	 * @return
	 */
	public int getRank()
	{
		Matrix xxM = new Matrix(getData());
		return xxM.rank();
	}

	/**
	 * Join two matrix according to qr.qy format (R function)
	 * 
	 * @param rr
	 * @param hh
	 * @return
	 */
	private static Mat triangularJoin(double[][] rr, double[][] hh)
	{
		int rows = hh.length;
		int cols = rr.length;
		double[][] triData = new double[rows][cols];
		for (int col = 0; col < cols; col++)
		{
			for (int row = 0; row < rows; row++)
			{
				if (col >= row)
				{
					triData[row][col] = rr[row][col];
				}
				else
				{
					triData[row][col] = hh[row][col];
				}
			}
		}
		Mat tri = new Mat(triData);
		return tri;
	}

	/**
	 * Return a matrix [rows,cols] full of value n
	 * 
	 * @param rows
	 * @param cols
	 * @param n
	 * @return
	 */
	public static Mat fillNmatix(int rows, int cols, double n)
	{
		double[][] data = new double[rows][cols];
		for (int i = 0; i < rows; i++)
		{
			double[] row = data[i];
			for (int j = 0; j < cols; j++)
				row[j] = n;
		}

		return new Mat("fillNmatrix(" + rows + "," + cols + "," + n + ")", data);
	}

	/**
	 * Return a matrix with only diagonal values and Zero elsewhere
	 */
	public Mat getDiagonal()
	{
		Mat diagonal = Mat.fillNmatix(rowDimension, colDimension, 0.0);
		for (int i = 0; i < colDimension; i++)
			diagonal.set(i, i, get(i, i));
		return diagonal;

	}

	/**
	 * Multiply for Matrix m, it is equivalent of operator <b>%*%</b> in R
	 * 
	 * @param m
	 * @return
	 */
	public Mat scalarMultiply(Mat m)
	{
		// cols A == row B ???;
		double val[][] = new double[rowDimension][m.colDimension];
		for (int i = 0; i < m.colDimension; i++)
		{
			double[] row = m.getColumn(i);
			for (int j = 0; j < rowDimension; j++)
			{
				double[] col = getData()[j];
				val[j][i] = scalarProduct(row, col);
			}
		}
		return new Mat(val);
	}

	/**
	 * Scalar product between two vectors a and b
	 * 
	 * @param a
	 * @param b
	 * @return
	 */
	private double scalarProduct(double[] a, double[] b)
	{
		double res = 0;
		for (int i = 0; i < a.length; i++)
		{
			res += a[i] * b[i];
		}
		return res;
	}

	/**
	 * Return in a matrix the sum of every column powered <b>exp</b>
	 * 
	 * @param exp
	 * @return
	 */
	public Mat sumColumnWisePower(int exp)
	{
		double[][] sum = new double[1][colDimension];

		for (int i = 0; i < rowDimension; i++)
		{
			double[] a = data[i];
			double[] c = sum[0];
			for (int j = 0; j < colDimension; j++)
				c[j] += Math.pow(a[j], exp);
		}

		return new Mat("sumColumnWisePower (" + name + ")^" + exp, sum);
	}

	public Mat toSquareForm()
	{
		int[] rows = new int[colDimension];
		for (int i = 0; i < colDimension; i++)
			rows[i] = i;
		return selectRows(rows, 0);
	}

	/**
	 * Read from a CSV files
	 * @param name
	 * @param fileName
	 * @param rows
	 * @param cols
	 * @param delimiter
	 * @return
	 * @throws IOException
	 * @throws InterruptedException
	 */
	public static Mat readFromCsvArray(String name, String fileName, int rows, int cols, String delimiter) throws IOException, InterruptedException
	{
		CsvMatrixReader reader = new CsvMatrixReader(name, fileName, rows, cols, delimiter);
		reader.process();
		Mat m = reader.getMatrix();

		if (verbose)
			m.print();

		return m;
	}

	/**
	 * return a matrix with incremental sum for COLS, columns
	 * 
	 * @param m
	 * @return
	 */
	public Mat cumsum()
	{
		return cumsum(COLS);
	}

	/**
	 * return a matrix with incremental sum form COLS, columns or ROWS, rows
	 * 
	 * @param m
	 * @param mode
	 * @return
	 */
	public Mat cumsum(int mode)
	{
		if (mode == ROWS)
			return cumsumRows();
		if (mode == COLS)
			return cumsumColums();
		return null;
	}

	private Mat cumsumColums()
	{
		double[][] r = new double[rowDimension][colDimension];
		double[][] m = getData();
		for (int row = 0; row < rowDimension; row++)
		{
			double[] m0 = m[row];
			double[] r0 = r[row];
			r0[0] = m0[0];
			for (int col = 1; col < colDimension; col++)
			{
				r0[col] = r0[col - 1] + m0[col];
			}
		}
		return new Mat("cumsum Colums(" + getName() + ")", r);
	}

	private Mat cumsumRows()
	{
		double[][] r = new double[rowDimension][colDimension];
		double[][] m = getData();
		r[0] = m[0];
		for (int row = 1; row < rowDimension; row++)
		{
			double[] m1 = m[row];
			double[] rr0 = r[row - 1];
			double[] rr1 = r[row];
			for (int col = 0; col < colDimension; col++)
			{
				rr1[col] = rr0[col] + m1[col];
			}
		}
		return new Mat("cumsum Rows (" + getName() + ")", r);
	}

	/**
	 * return the NOT (1 instead 0 and vice versa) matrix from a 0/1 matrix
	 * 
	 * @return
	 */
	public Mat not()
	{
		double[][] r = createResult();
		for (int row = 0; row < rowDimension; row++)
		{
			double[] a = r[row];
			for (int col = 0; col < colDimension; col++)
			{
				a[col] = (get(row, col) == 0.0 ? 1.0 : 0.0);
			}
		}
		return innerOrNew("not (" + getName() + ")", r);

	}

	/**
	 * The difference between two columns
	 * 
	 * @return
	 */
	public Mat diff()
	{
		int rows = rowDimension;
		int cols = colDimension - 1;

		Mat m = new Mat("diff", rows, cols);
		double[][] that = m.getData();

		for (int i = 0; i < rows; i++)
		{
			double[] a = data[i];
			double[] b = that[i];
			for (int j = 0; j < b.length; j++)
			{
				b[j] = a[j + 1] - a[j];
			}
		}

		return m;
	}

	/**
	 * return a mat sorted by ROWS or COLS
	 * 
	 * @param mode
	 * @return
	 */
	public Mat sort(int mode)
	{
		if (mode == Mat.ROWS)
			return sortRows();
		if (mode == COLS)
			return sortCols();
		return null;
	}

	private Mat sortCols()
	{
		double[][] b = new double[colDimension][rowDimension];
		for (int col = 0; col < colDimension; col++)
		{
			double[] a = getColumn(col);
			Arrays.sort(a);
			b[col] = a;
		}
		Mat res = new Mat("Sort COLS(" + getName() + ")", b);
		res = res.transpose();
		return res;
	}

	private Mat sortRows()
	{
		double[][] b = new double[rowDimension][colDimension];
		for (int row = 0; row < rowDimension; row++)
		{
			double[] a = getRow(row);
			Arrays.sort(a);
			b[row] = a;
		}
		Mat res = new Mat("Sort ROWS(" + getName() + ")", b);
		return res;
	}
	
	
	public void sortByColumn(Comparator<double[]> comparator)
	{
		Arrays.sort(data, comparator);		
	}
	
	
	public void sortByColumn(final int column)
	{
		sortByColumn(new Comparator<double[]>() {

			@Override
			public int compare(double[] o1, double[] o2)
			{
				double d1 = o1[column];
				double d2 = o2[column];
				
				return (int)Math.signum(d2 - d1);
			}		
		});
	}


	/**
	 * return the Minor matrix for Row rr and Column cc
	 * return the matrix that NOT contain the Row rr and the Rolumn cc
	 * @param rr
	 * @param cc
	 * @return
	 */
	public Mat minorMatrix(int rr, int cc)
	{
		double[][] res = new double[rowDimension - 1][colDimension - 1];
		
		int rowIndex = 0;
		for (int r = 0; r < rowDimension; r++)
		{
			if (r == rr)
				continue;
			
			double[] dataRow = data[r];
			double[] resRow = res[rowIndex];
			
			int colIndex = 0;
			for (int c = 0; c < colDimension; c++)
			{
				if (c == cc)
					continue;
				
				resRow[colIndex++] = dataRow[c];
			}
			rowIndex++;
		}
		return new Mat("Minor(" + name + ")", res);
	}
	
	/**
	 * return the adjugate of a matrix
	 * @return
	 * @throws Exception
	 */
	public Mat adjugate() throws Exception
	{
		if (rowDimension != colDimension)
			throw new MatException("Matrix Not Square");

		double[][] res = new double[rowDimension][colDimension];
				
		MatRowWorker job = new MatRowWorker("", res, res, null) {
			@Override
			public int proposeNumberOfThreads()
			{
				return 100;
			}

			public void run()
			{
				for (int i = iStart; i < iEnd; i++)
				{
					double[] row = result[i];
					for (int c = 0; c < colDimension; c++)
					{
						double exponent = ((i + c) % 2 == 0 ? 1.0 : -1.0);
//						Mat b = minorMatrix(i, c);
						
						int[] ir = drop(mkidx(rowDimension),i);
						int[] ic = drop(mkidx(colDimension),c);
						
						row[c] = exponent * _det(ir, ic);
					}
				}
			}
		};
		Mat m  =  executeMatRowWorkers(job);
		m = m.transpose();
		m.setName(concat("Adjugate(", name, ")"));
		
		return m;
		
		
//		return new Mat(concat("Adjugate(", name, ")"), res).transpose();
	}

	// too slow, goes by n!
//	public Mat inverseMatrix() throws Exception
//	{
//		double det = det();
//		Mat inv = adjugateMatrix();
//		inv = inv.inner().scalarMultiply(1 / det);
//		inv.setName(concat("inverse(",inv.getName(),")"));
//
//		return inv;
//	}
	
	/**
	 * inverses matrix using Gauss Jordan algorithm
	 */
	public Mat inverse() throws MatException
	{
		if (rowDimension != colDimension)
			throw new MatException("Matrix Not Square");

		Mat matA = this.copy("A");
		Mat matB = Mat.diag(concat("inverse(",name,")"), rowDimension, 1); 
		
		return gaussJordan(matA, matB);
	}

		
	private Mat gaussJordan(Mat matA, Mat matB) throws MatException
	{	
		double[][] a = matA.getData(); 
		double[][] b = matB.getData(); 
		int n = b.length; 
		int colA = matA.colDimension;
		int colB = matB.colDimension;
		
		for (int k = 0; k < n; k++)
		{
			pivot(a,k,b);
			double[] ak = a[k];
			double[] bk = b[k];
			double d = ak[k];
			if (d == 0.)
				throw new MatException("Singular matrix, row/col=" + k);
				
			for (int j = 0; j < colA; j++)
				ak[j] /= d;
			
			for (int j = 0; j < colB; j++)
				bk[j] /= d;

			
			for (int i = 0; i < n; i++)
			{
				if (i == k)
					continue;
				
				double[] ai = a[i];
				double[] bi = b[i];
				d = ai[k];
				
				for (int j = 0; j < colA; j++)
					ai[j] -= d*ak[j];
				
				for (int j = 0; j < colB; j++)
					bi[j] -= d*bk[j];
			}
		}
		
		return matB;
	}
	
	/**
	 * exchange row k with row whose row(i > k,k) is max. Used
	 * with gauss jordan to avoid diagonal element being 0 which
	 * would cause a "divide by zero" exception.
	 * 
	 * @param a	check for max at col x
	 * @param k
	 * @param b	swap a AND b
	 */
	private void pivot(double[][] a, int k, double[][] b)
	{
		int kMax = k;
		double max = Math.abs(a[k][k]);
		for (int i = k + 1; i < a.length; i++)
		{
			double d = a[i][k];
			if (Math.abs(d) > max)
			{
				kMax = i;
				max = d;
			}
		}
		
		if (kMax != k)
		{
			// swap
			double[] tmp = a[k];
			a[k] = a[kMax];
			a[kMax] = tmp;	
			
			tmp = b[k];
			b[k] = b[kMax];
			b[kMax] = tmp;			
		}
	}

	
	/**
	 * Return the determinant of a matrix
	 * @return
	 * @throws Exception
	 */
	public double det() throws Exception
	{
		if (rowDimension != colDimension)
			throw new MatException("Matrix Not Square");
		
		return _det(mkidx(rowDimension), mkidx(colDimension));
	}	
	
	/**
	 * Return a new matrix without the rows specified into rows
	 * @param rows
	 * @return
	 */
	public Mat deteteRows(Integer... rows){
		ArrayList<Integer> rowsToDelete = new ArrayList<Integer>(Arrays.asList(rows));
		double[][] res = new double[rowDimension-rows.length][colDimension];
		int index =0;
		for(int row=0;row<rowDimension;row++){
			if(rowsToDelete.contains(row))
				continue;
				res[index++]=data[row];
		}
		return new Mat(name,res);		
	}
	
	/**
	 * Return a new matrix without the columns specified into cols
	 * @param rows
	 * @return
	 */
	public Mat deleteCols (Integer... cols){
		Mat a =this.transpose();
		a = a.deteteRows(cols);
		return a.transpose();
	}
	

	public Mat skipNan(int mode) {
		if (mode == Mat.ROWS)
			return skipNanRows();
		if (mode == COLS)
			return skipNanCols();
		return null;
	}
	
	
	
	
	/**
	 * Delete the rows that contains any NaN
	 * @return
	 */
	private Mat skipNanRows()
	{
		ArrayList<Integer> rowsWithNaN = new ArrayList<Integer>();
		for(int r =0; r<rowDimension;r++){
			double[]row= data[r];
			for(int c =0 ; c<colDimension;c++){
				if(new Double(row[c]).isNaN())
					rowsWithNaN.add(r);
			}
		}
		if(rowsWithNaN.size()==0){
			return this;
		}
		else
		{
			return this.deteteRows(rowsWithNaN.toArray(new Integer[rowsWithNaN.size()]));
		}
	}
	
	/**
	 * Delete the cols taht contains any NaN
	 * @return
	 */
	private Mat skipNanCols()
	{
		Mat a = transpose();
		a=a.skipNanRows();
		return a.transpose();
	}

	private double _det(int[] ir, int[] ic)
	{	
		double res = 0;
		if (ir.length == 2)
		{
			double[] r0 = data[ir[0]];
			double[] r1 = data[ir[1]];
			int c0 = ic[0];
			int c1 = ic[1];
			res = r0[c0] * r1[c1] - r0[c1] * r1[c0];
		}
		else if (ir.length == 3)
		{
			double[] r0 = data[ir[0]];
			double[] r1 = data[ir[1]];
			double[] r2 = data[ir[2]];
			int c0 = ic[0];
			int c1 = ic[1];
			int c2 = ic[2];

			res = r0[c0] * r1[c1] * r2[c2] + r0[c1] * r1[c2] * r2[c0]
					+ r0[c2] * r1[c0] * r2[c1] -

					r0[c2] * r1[c1] * r2[c0] - r0[c1] * r1[c0]
					* r2[c2] - r0[c0] * r1[c2] * r2[c1];
		}
		else
		{
			double[] row1 = data[ir[1]];
			int[] row1Dropped = drop(ir,1);
			boolean add = true;
			for (int c = 0; c < ic.length; c++)
			{
				double cell = row1[ic[c]];
				if (cell == 0.)
					continue;
				add = !add;
				
				// minorMatrix(1, c);
				double d = _det(row1Dropped, drop(ic,c));
				if (add)
					res = res + cell * d;
				else
					res = res - cell * d;
			}
		}
		
		return res;
	}
	
	private static int[] mkidx(int size)
	{		
		int[] a = new int[size];
		
		for (int i = 0; i < a.length; i++)
			a[i] = i;
		
		return a;
	}
	
	private static int[] drop(int[] a, int idx)
	{
		int[] res = new int[a.length - 1];
		
		int j = 0;
		for (int i = 0; i < a.length; i++)
		{
			if (i != idx)
				res[j++] = a[i];
		}
		
		return res;
	}

	public Mat divideColumn(double[] b)
	{
		int rows = rowDimension;
		int cols = colDimension;
		double[][] result = createResult();
		for(int row =0;row<rows;row++){
			double[] a = data[row];
			double[] c = result[row];
			for (int col= 0; col < cols;col++)
					c[col] = a[col] / b[col];
		}
		return innerOrNew(concat("(", this.name, ".divideColumn", ")"), result);
	}
	
	
	
	
	

}


